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ABSTRACT 

Aims. 3D tomography of the interstellar dust and gas may be useful in many respects, from the physical and chemical evolution of 
the ISM itself to foreground decontamination of the CMB, or various studies of the environments of specific objects. However, while 
spectral data cubes of the galactic emission become increasingly precise, the information on the distance to the emitting regions has 
not progressed as well and relies essentially on the galactic rotation curve. Our goal here is to bring more precise information on the 
distance to nearby interstellar dust and gas clouds within 250 pc. 

Methods. We apply the best available calibration methods to a carefully screened set of stellar Stromgren photometry data for targets 
possessing a Hipparcos parallax and spectral type classification. We combine the derived interstellar extinctions and the parallax 
distances for about 6,000 stars to build a 3D tomography of the local dust. We use an inversion method based on a regularized 
Bayesian approach and a least squares criterion, optimized for this specific data set. We apply the same inversion technique to a 
totally independent set of neutral sodium absorption data available for about 1700 target stars. 

Results. We obtain 3D maps of the opacity and the distance to the main dust-bearing clouds with 250 pc and identify in those maps 
well-known dark clouds and high galactic more diffuse entities. 

We calculate the integrated extinction between the Sun and the cube boundary and compare with the total galactic extinction derived 
from infrared 2D maps. The two quantities reach similar values at high latitudes, as expected if the local dust content is satisfyingly 
reproduced and the dust is closer than 250 pc. Those maps show a larger high latitude dust opacity in the North compared to the South, 
reinforcing earlier evidences. Interestingly the gas maps do not show the same asymmetry, suggesting a polar asymmetry of the dust 
to gas ratio at small distances. 

We compare the opacity distribution with the 3D distribution of interstellar neutral sodium resulting from the inversion of sodium 
columns. We discuss the similarities and discrepancies and the influence of data set limitations. Finally we discuss the potential 
improvements of those 3D maps. 
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1. Introduction 

The 3D representation of the local ISM is an important chal- 
lenge in various perspectives. In the most direct way it allows a 
better modelling of the ISM. The ability to localise masses of 
gas and dust, combined with information on the ionisation state, 
the temperature and the kinematics helps to understand the 
chemical and physical evolution of the different gas components 
in relation with their history and the surrounding ionizing 
radiation from hot stars, as well as the dust-gas relationships 
and the dust processing (Cox, 2005; Draine, 2003). 
In a more indirect way, the 3D ISM structure becomes increas- 
ingly useful as an ingredient to the foreground contamination 
of the CMB because the local distribution of the gas influences 
the interstellar radiation field, which in turn influences the 
dust temperature and infrared-millimetric emission. The 3D 
distribution of the local ISM is also governing a significant 
fraction of the diffuse gamma-ray background and is one of the 
ingredients of the models developed in the context of the Fermi 
satellite. 

In a very different way, the knowledge of the 3D properties of 



the local ISM may be useful in various other studies since it 
provides the environmental context of astrophysical objects, al- 
lows a correction of the extinction according to the distance and 
the direction, allows to disentangle spectral lines in absorption 
related to the object itself or to the foreground, or provides a 
tool for the estimate of the propagation of energetic particles. 
The list of topics is too long to be detailed here. 

Recent work on the 3D global structure and kinematics of 
the local ISM or on specific nearby ISM structures has been 
done using absorption lines in the UV (e.g. Lallement et al, 
1995, Redfield and Linsky, 2008, Welsh and Lallement, 2005) 
and the strong Nal and Call interstellar in the visible (e.g. 
Genova and Beckman, 2003, Snow et al, 2008). A number of 
works have combined emission data and neutral sodium ab- 
sorption to attribute a distance to specific features seen in ra- 
dio or IR (e.g. Lilienthal et al, 1992, Corradi et al, 2004, Meyer 
et al, 2006), making use essentially of Hipparcos parallaxes 
(Ferryman et al.(I997) i. Neutral species like Nal trace low tem- 
perature regions, i.e. dense atomic gas clouds and molecular 
clouds, while ionised calcium is a tracer of both dense regions 
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(at the exception of the fully neutral cores) and warm diffuse 
gas, partially ionised (see Welty et al, 1996 for more details). 
Note that the known existence of different typical sizes linked 
to the different types of gas, from parsecs for molecular clouds 
to tens of parsecs for atomic and diffuse clouds, has influenced 
our choice of correlation lengths for the inversion described in 
section 4. In parallel, Stromgren photometry has been used in 
combination with parallaxes for distance estimates of a number 
of features (e.g. Reis et Corradi , 2008, Knude and Hog 1998 , 
1999) 

It is known for a long time from both absorption data 
and diffuse soft X-ray background detection that the Sun is 
embedded in a relatively empty cavity, the Local Bubble (LB), 
mostly devoid of dense clouds (Frisch and York, 1983, Sfeir 
et al, 1999). Within this local cavity reside essentially tenuous 
diffuse clouds such as the clouds that make part of the so-called 
Local Fluff within 20 pc from the Sun. The interpretation of the 
soft X-ray background data suggests that the cavity is filled by 
one million K gas (Snowden et al, 1990). The cavity is delimited 
by dense clouds whose distances range from 40 to 180 parsecs 
and neighboring cavities are distributed all around. Such a 
"cellular" structure of hot bubbles separated by high density and 
relatively compact regions, including molecular dust clouds is 
totally expected from the permanent ISM recycling that shapes 
it: stellar winds and supernova explosions inflate hot "bubbles" 
within the dense gas and maintain these high contrasts (Cox and 
Reynolds, 1987, De AviUez and Breitschwerdt, 2004). There is 
very probably a strong link between the strong winds having 
inflated the cavity and the surrounding so-called Gould belt, 
an expanding elliptical ring of star-forming regions (see Perrot 
and Grenier, 2003). Several questions remain however about 
the local ISM and Local Bubble. The embedded clouds are 
ionized by the early-type stars, but their measured ionisation 
state, especially the ionization of helium requires an additional 
harder radiation supposed to be generated within semi-hot 
conductive interfaces between clouds and hot gas (e.g. Slavin 
1989). The distribution of the corresponding semi-hot or hot 
gas ions however is not in fully satisfying agreement with 
this scenario (Welsh and Lallement, 2005, Savage and Lehner, 
2006, Knauth et al, 2003). On the other hand the more recent 
discovery of solar wind charge-transfer X-ray emission (Cox, 
1998, Cravens, 2000) has also complicated the picture since 
this diffuse emission, which mimics one million K thermal 
emission, is as intense as the observed background in a number 
of directions (e.g. Koutroumpa et al, 2009). Decreasing the 
hot gas pressure significantly would have the positive effect of 
suppressing the strong difference from the pressure measured 
within the diffuse clouds and the local cloud (Jenkins, 2009). 
From 3D gas and dust mapping one can expect significant 
progresses on these topics, the "Local Bubble" wall location 
and the X-ray emission source regions, some hints about 
the local bubble formation and the link between the local 
bubble and surrounding bubbles and more generally a better 
understanding of the multi-phase structure and the interaction 
between the different phases of the ISM. One can also expect 
the determination of the regions of star formation and estimates 
of the cloud masses. Last but not least, it should help to shed 
light on the surprisingly high deuterium abundance in the local 
ISM (Linsky et al, 2006), now largely recognised to be related 
in a very large extent to dust evaporation from grains (Draine, 
2004). The unambiguous correlation between refractory metals 
and deuterium abundances in the gas (Lallement et al, 2008) as 
well as the peculiar distribution of the deuterium enhancements 
with respect to the Gould belt (Lallement, 2009) suggest a 



scenario for the formation of the belt with a large role of dust 
and dust-gas processes. 

The first computed tomography of the local ISM from 
[Vergely et al.(200r)| used neutral sodium absorption data com- 
bined with Hipparcos stellar parallaxes (Ferryman et al, 1997) 
and a tomographic method based on the work of Tarantola & 
Valette (1982a), specifically applied for the first time to extinc- 
tion by the ISM by |Chen et al.(19 98) This type of tomography 
is a new manner to represent the ISM, considering it as a contin- 
uous medium without a priori cloud discretization and consid- 
ering the clouds and the voids as fluctuations around a median 
density, the "prior" value. Lallement et al.(200 3) presented ad- 
ditional sodium observations especially dedicated to the Local 
Bubble boundary determination, recorded during several observ- 
ing programs, and made use of about 1,000 carefully screened 
data for a column density mapping and a tomography based on 



the same algorithm as in Vergely et al.(2001)| An updated study 
using about 1700 targets as well as the first tomography using 
ionized calcium are presented in the recent study by Welsh et al 
(2009). Increasing the number of targets has allowed to reveal 
neighboring cavities and to gain in precision. 
Using a similar approach, we present here the 3-D density dis- 
tribution of the extinction in the vicinity of the Sun, using color 
excesses from Stromgren photometry and Hipparcos distances. 
Because extinction traces low temperature areas as does neutral 
sodium, a comparison of the derived extinction maps with the 
results of Welsh et al. 2009 is performed. 

Section 2 describes in detail the determination of the extinc- 
tion and the associated uncertainties. Section 3 describes the ba- 
sis of the inversion method and the algorithm used in our study as 
well as the underlying assumptions. Section 4. 1 compares the to- 
tal extinction towards the target stars with the total extinction de- 
duced from IR data by Schlegel et al (1998), and compares inte- 
grated extinctions and Nal columns. Section 4 shows the results 
of the inversion in several planes and the comparison with the in- 
version of the neutral sodium absorption. Section 5 presents the 
conclusions and discusses shortcomings and potential improve- 
ments of the inversion and the database. 

2. Photometric data and extraction of extinction 

2. 1 . Extinction representation 

In the standard model of stellar light absorption by intervening 
dust clouds, the intensity decreases exponentially along the line- 
of-sight as a function of the opacity : 

/(i) = /o(/l)exp(-T,0 (1) 

1(A), /()(/!) are the absorbed and unabsorbed intensities and is 
the optical depth, or opacity, depends on the wavelength A, the 
absorption coefficient k,i and the density of the dust p^,,.,,, and is 
the integral along the line-of-sight (LOS) up to the distance Ri 
of the star: 
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ki(r)pdua(r)r 



(2) 



The observed flux decreases with the dust absorption and 
with the distance. To quantify the dust absorption solely and re- 
move a dependence on the intrinsic stellar flux and the distance 
one uses the color excess Ei(Ai,A2) that is the logarithm of the 
intensity ratio between two different bands Ai and A2 : 
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EiiAi , iz) = C (k^, (r) - k^^ (r)) Pd„„(r)r 



(3) 
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Fig. 1. Histogram extinction for stars in the 70 pc sphere 



In the following sections we adopt the notation : 

Ei Ei(Ai,A2) and k := - k^. 



(4) 



Note that k,i(r) depends on the composition of the dust which 
can vary along the line of sight. However our purpose here is to 
estimate the opacity and not to extract the actual dust density 
by using an extinction coefficient. Thus, the variable of interest 
is the differential opacity, a combination of the dust properties 
(dust size distribution and extinction coefficient) and the dust 
density: 



Pop(r) := C (^i,(r) - k^^ir)) pdust(r) 
In this case 



Ei(AuA2)= Pop{r)r 
Jo 



(5) 



(6) 



In this paper we use the differential opacity between the 
Stromgren photometry B and y bands: 

Popir) := C (^b('") - K{r))pdust{r) (7) 
The quantity pop is expressed in color excess per pc. 



The tomography of the extinction assumes a continuous dis- 
tribution of the opacity in 3D and requires a 3D formulation of 
Eqn. |6] For each line-of-sight / among the n measured opacities 
constituting the data set, the /th constraint is expressed in the 
following way over the whole space V: 



£, = C Ki(x)k{x)pdust{x)y{x)^ Ux)pop{x)y{x) (8) 



= r u 

Jv 



where x is the position vector in galactic coordinates : x = 
(r,l,b) with y(x) = r^cos(b)Thl, and where the singular kernel 
Ki(x) is explicitly given by : 



Ki(r,l,b)^@iir)6ib-bi)6il-ld 



(9) 



in denoting by 0,(r) the function equal to l/r^cos{bi) inside the 
interval [0, /?,] and null outside, and by 6 the Diracs' distribution. 

The density of IS matter decreases strongly with the distance 
to the galactic plane. A convenient description of this behaviour 
is to suppose a priori an exponential decrease with height above 
the galactic plane (Chen et al., 1998) : 



Pref(.r,b) =poexp 



|r sin(/7)| 



(10) 



4 



J-L. Vergely, B. Valette, R. Lallement, and S. Raimond,: Spatial distribution of interstellar dust in the Sun vicinity 



We assigned a value of 200 pc to ho for both the opacity and Nal. 
Note that for dust and gas 200 pc is somewhat larger than what 
is usually used (100 pc for the dust, see, i.e., Chen et al., 1998, 
and 160 pc for Nal as deduced from a fit to the column-densities 
of our Nal database), but we have done it intentionally in order 
to avoid the disappearance of the tenuous high latitude structures 
along the inversion when the "a priori" solution is too small. 
In this study, we propose to find an extinction model close to 
the reference model given in Eqn. 10 that explains observed star 
extinctions. 

Since the density is a positive quantity, it is convenient to make 
the following change of variable : 



p(x) = pref exp (a(x)) 



(11) 



where the log opacity a is the new unknown parameter, which 
is assumed to follow a gaussian centered law, and which traces 
the difference in opacity from the reference model. In the 
following numerical applications, we have taken an average 
E(b-y) color excess po = 0.0004 per pc. This choice, which 
corresponds to about twice the average measured value, is 
motivated by the following reasons: i) we aim at a robust result 
on the very low opacity within the Local Bubble. Obtaining 
a strong decrease from the "a priori" model demonstrates the 
observational constraints on the cavity, ii) we aim at showing 
here that the inversion algorithm does not depend on the " a 
priori" opacity, iii) in areas with high dust concentration there 
may be an observational bias towards weaker opacities due 
to the limited sampling (lines of sight towards the targets do 
not cross the densest cores). Starting with larger opacities may 
partially compensate for this effect. 



2.2. Extinction determination from Stromgren ptiotometry 

For the derivation of color excesses we use the Stromgren color 
indices from the Hauck-Mermilliod Catalogue ( 1990). These au- 
thors have produced an extensive compilation of photometric 
data, and in case of multiple measurements they have computed 
the averaged values, taking into account the individual errors. 
The catalog lists 66,000 objects. We have selected from the cat- 
alog those stars possessing H/S determinations and Hipparcos 
parallaxes higher than 3.34 mas (about 16,000 stars). We have 
also very carefully screened the data using the Tycho catalog 
and removed all stars suspected or observed to be part of a mul- 
tiple system, all variable stars and all stars suspected to be sur- 
rounded by a shell. Using the Tycho spectral type classification 
all luminosity class I and II have finally also been removed. The 
remaining data set (about 6400 targets) has then been entered in 
the calibration process. 

Using calibration relations and spectral type classification, 
color excesses E{b - y) were derived for four groups of stars. 
Table 1 gives for the spectral type and the range of /3 index for 
the different calibration relations used in this study. A detailed 
description of the calibrations is given in [Vergely et al.(2001) 



Because the tomographic technic is very sensitive to the pres- 
ence of outliers, a specific filter has been applied on color ex- 
cesses from group 1. As mentioned by Crawford ( 1979), target 
stars with high cl values (evolved stars with high luminosity) 
should be removed. This is due to the presence of a turn off in 
the (cl,(b-y)) color diagram which generates an ambiguity in 
the location of the zero age main sequence. Indeed, a careful in- 
spection of stars with cl > 0.9 shows clearly the presence of 
underestimated E{b - y) color excesses and we have removed 
the corresponding targets. Note that the threshold recommended 
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Fig. 2. Fit between modeled data and observed one. 



0.35 



from Crawford (1 19791 1 is 0.58 instead of 0.9, which may ques- 
tion our use of the class III stars. However only 5 stars (from 
6,400) do not fulfil this limitation. We have checked that those 
stars, which are found to have a very small color excess, show 
no disagreement with the surrounding targets. 

Assuming that the expected E(b-y) color excess is negligible 
in the 70 pc sphere, it is possible to determine, for each calibra- 
tion group, the intrinsic dispersion and the bias of the observed 
color excesses. The figure [T] shows that only slight biases affect 
the different group. In order to avoid a propagation of such a bias 
in the tomography output, this bias has been corrected. The ta- 
ble 2 gives the bias and the standard deviation of the color excess 
from star belonging to the 70 pc sphere. The standard deviation 
comes from intrinsic dispersion in color-magnitude diagram and 
photometry error. It gives us the possibility to estimate the color 
excess error to be applied during the inversion. 
We have checked our results on the extinction based on the 1997 
Hauck Mermilliod Stromgren data base by comparing them to 
previous results by Philip and Egret (1980). These authors have 
used similar calibration relations for group 1, 3 and 4 and the 
earlier Hauck-Mermilliod (1980) catalogue. The cross correla- 
tion between our sample and the Philip and Egret sample allows 
to extract 2456 common stars. The standard deviation between 
the differential opacities E(b-y) from Philip and Egret and those 
from the present study reaches 0.009 mag and the mean differ- 
ence is -0.003 mag. Most of those differences are due to differ- 
ences in the Stromgren photometry colors which are more recent 
in the present study and in differences in calibration relations (for 
intermediate group, we use Hilditch et al.(1983) calibration). 

2.3. Target distances and distribution 

The distance to the stars affected by extinction in the solar neigh- 
borhood is given by the Hipparcos parallaxes. In this study, we 
use a new set of Hipparcos parallaxes consistently derived by 
IVa n Leeuwen (2007) | This new Hipparcos processing provides 
parallaxes which are slightly different but of importance for our 
study is the resulting accuracy, claimed to be a factor of 2.2 bet- 
ter than the initial one (from Ferryman et al.(1997)") . Because the 
distance accuracy reaches about 100 pc at 400 pc, we have kept 
stars closer than 300 pc, hence our limit on the parallax quoted 
above. 

The distribution of retained targets is far from regularly dis- 
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group 


spectral type 


luminosity class 




number of stars 


calibration 


1 


BO to AO 


III to V 


2.59 </?< 2.88 


993 


Crawford (1978) 


2 


AO to A3 


III to V 


2.87 </?< 2.93 


626 


Hilditch et al.(1983) 


3 


A3 to FO 


nito V 


2.12 < 13 < 2.88 


1517 


Crawford (1979) 


4 


FO to G2 


m to V 


2.56 </?< 2.72 


3857 


Crawford (1975) 


TOTAL 








6993 





North 



South 




0.15 0.2 0.25 
E(b-y) (this study) 



Fig. 3. Comparison between Schlegel et al. extinction and ex- 
tinction obtained in this study ; north galactic pole in Lambert 
projection. 



tributed in 3D. On the contrary, there are areas of high concentra- 
tion, and generally there is a strong drop of target density beyond 
about 100 parsecs. The distribution of targets is illustrated in the 
figures presented in section 5. 

2.4. Error estimates 

If observed data are biased or if their errors have been underes- 
timated, we cannot expect to obtain significant results because 
the inversion technique tries to reproduce an image of the initial 
data. Data errors arise both from extinction and distance errors. 
In our model, we do not consider these two errors as independent 
but we propagate the distance error on that of the extinction. The 
distance error o-j is related to the error of the parallaxes, cr„ in 
the following way : 

1^1 = 1^1 
a n 

The extinction error crg^ resulting from the distance error is ap- 
proximately : 

ctej = E(b - y) ^ 
a 

if we assume that the opacity is constant along the line of sight. 
So, if we consider a residual error of cTcai for the extinction cali- 
bration, the total error, ctej-, will be in the gaussian case: 



cal 



Table |2] shows the calibration error deduced from the standard 
deviation of the stars belonging to the 70 pc sphere. 




0.15 0.2 0.25 
E(b-y) (this study) 



Fig. 4. Same as Fig.^for the south galactic pole. 
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Fig. 5. Sodium column density versus E(b-y) for common line 
of sight. 



3. Tomography of the opacity 

3.1. The inverse method 

The determination of the 3-D scalar field from discrete obser- 
vational data is clearly under constrained because the sample 
of stars is limited to a finite number of lines of sight rela- 
tively spread-out. The problem can be regularized by impos- 
ing that the true solution is smooth in some suitably quantified 
sense (Twomey 1977; Craig & Brown 1986). Thus one reduces 
the number of free parameters, and the problem becomes well- 
posed. 
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Table 2. Estimation of errors and biases on extinction according to the group of calibration, for stars in the 70pc sphere 



group 


Ij^rror 


Bias 


1 


0.0044 


-0.0024 


2 


0.0145 


-0.0011 


3 


0.0088 


-0.0009 


4 


0.0117 


0.0040 



Here, we follow the approach of Tarantola & Valette 
(1982a,b). In this stochastic approach the prior information is 
represented by probabihty measure for both data and model pa- 
rameters. Thus true data d and model m are considered as ran- 
dom vectors to take account of the errors in the data, in the one 
hand, and of the "a priori" information on the model, in the other 
hand. The inversion process consists in imposing that the data 
are in fact linked to the model through a theoretical mapping 
g : d = g(m). This is achieved by considering, as posterior in- 
formation on the model, the conditionnal probabihty law for m 
given that the a priori random vector d - g{m) is in fact null. In 
the present case of extinction, m is the opacity (or the log opac- 
ity) function over the space, and g is given by (eqn. 8). When 
both data and model are assumed to be Gaussian random vec- 
tors, with centre and covariance operator respectively denoted 
by dprior, Cd, and Mprwr, Cm, and when g is weakly non-linear, 
Tarantola and Valette (1982b) showed that the posterior measure 
for m can be well approximated by the Gaussian measure, the 
centre m of which mjikes minimum E{m): 

E(m) = \\C-l'Hm - mp,,or)f + \\C-/'\g(m) - dprior)f (12) 

and the covariance operator of which is expressed as : 

C = (C-J +GmC-/GJ-' (13) 

where G*^ denotes the adjoint - with respect to the usual scalar 
product of the model space and the data space - of the derivative 
G„i of g at model m. The model m verifies the following equation 
(Tarantola and Valette, 1982b) : 

^ - fnprior = CmG*^{Cd + GmCi„G*f^) ^ 

(dprior - g{m) + G^(m - nipnor)) (14) 

which can be solved iteratively (Tarantola and Valette, 1982b) 
following a fixed point algorithm : 

nik+l - mprior = CmG\{Cd + GkCmGl)'^ 

(dprior - g{mk) + Glcimk - ntprior)) (15) 

where Gk denotes the derivative of g at model m^. It is important 
to note that even though the model (i.e. the log-opacity a) is a 
function of the space, an iteration only requires to solve a linear 
system with rank equal to the number of data. From a practi- 
cal point of view, the covariance operator of the log-opacity is 
defined as / ^ C^f with : 

Cmfix) = (r{x) cr(x') ,J,{x, x') f{x') Y(x') (16) 

where >J/(x, x!) is the correlation of the log-opacity between point 

X and x', and where cr^(x) is its variance at point x. Hence, it 
follows that the algorithm (15) takes the following simple form : 

ak+\{x) = (t{x) y 6' I p{Xi) (T{Xi) ij/{x, Xi) r,. (17) 



where n is the number of data, and where the n- vector 6 is solu- 
tion of the following linear system : 

S ij 6j = Ei+ [' pkixi) (akixi) - 1) r,. (18) 
Jo 

with 

S ij = {Cd)ii + a{Xi)pk{Xi) a{Xj)pk{Xj) il/{Xi, Xj) r,r/19) 

3.2. Resolving power and averaging index 

The concept of resolving kernel, introduced by Backus & Gilbert 
(1970), is an effective tool to appreciate the efficiency of inver- 
sion. Taking that d = g(p) into accoimt, a first order expansion 
of g around m - m yields 

dprior - g(m) + Ga(w - ntprior) = dprior -d + Gfh{m - nipriorX^O) 

Substituting this expression into (eqn 14) gives : 

m — Tn prior — Rm^ni — III prior) ~ Efii{d — dprior) (21) 

where R,;, is the resolution operator : 

Rm = CmG'fhiCd + GfiiC„G*^~^Gfh (22) 
and 

Em — CrriG-^{Cd + GfnCmGj^) , — E^Gm (23) 

Expression (21) shows that, correct to first order, the model m - 
Mprior is the sum of a term which corresponds to the difference 
between the true model and the prior m - niprior, filtered through 
the resolution operator, and of a noise term coming from the 
discrepancy between true and prior data. This expression can be 
rewritten as: 

(m - mprior){x) - I h{x, x') {m - mprior){x') V(x') -I- noise (24) 

Jv 

where h{x, x') is the kernel of the resolution operator, which is 
called the resolving, or the averaging kernel. Clearly, the value of 
m- niprior at point x is an average of the values of the true model, 
weighted by the kernel, all around, and the resolution would be 
perfect if this kernel was equal to the Dirac's distribution 6{x - 
xf). Given a point x, h(x, x') is generally a peaked function of x', 
centered around x' = x, with few negative lobes, and there are 
many ways to define its width (Rodgers, 2008, p. 54). 

Here, we will focus on another important and useful con- 
cept, closely related to the averaging kernel. Let (m - niprior) be 
the mean value of the true model within the width of the kemel 
around x. Equation (24) reads : 

(m - mprior)(x) = I{x) (m - ntprior) + noise (25) 
where : 

/(x):= r hix,x)Y(x') (26) 
Jv 
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is called the averaging index. In those areas where the data pro- 
vides poor information, I(x) takes very low values, and hence, 
the resulting model remains close to the prior one, independantly 
of the value of the true model in the neighbourhood. In contrast, 
where the index is close to one, the resulting model corresponds 
to the mean value of the true model within the width of the av- 
eraging kernel. Thus, displaying the area where the index is low 
is a way to indicate the limit of the zone where the inference 
becomes poor 

3.3. Prior information on tlie opacity distribution. 

In this method the model parameters may include finite dimen- 
sional vectors as well as functions. For finite dimensional model 
spaces. Cm is simply a matrix. In the case in which we do not 
have any prior information on parameters, the parameters must 
be a priori decorrelated with large variances. In the case of func- 
tions, we must define the correlation kernel t//(x, x') of the co- 
variance operator (16). Simple kernels are Gaussian kernels 



i//(x,x') - exp 



Wx-x'f 
or exponential kernel : 



i//(x, x') - exp 



^0 



(27) 



(28) 



where is a smoothing length that controls the spatial variabil- 
ity of the random function. 

In this study, we first used a Gaussian kernel that yielded a poor 
data fit. Since interstellar medium clearly shows different char- 
acteristic lengths, we found more judicious to use multiscale ker- 
nels (Serban and Jacobsen (2001)). Actually, the introduction of 
two scales, the first one characterizing the diffuse matter, and 
the second one, the compact clouds or group of clouds, was suf- 
ficient to obtain a good adjustment of data, in a second step. The 
CO variance kernel that we have used is the following: 



cr(x) cr(x') ipix, x) — cro(jic) croix') exp - 



cri(x) cri(x') exp 



\\x - X 



4i 



(29) 



The characteristic lengths and ^\ were taken respectively 
equal to 15 and 5 pc. The levels of fluctuation ctq and ct] were 
taken equal to 1.5 and 0.8 respectively. The choice of this dou- 
ble scale is important: it corresponds to the two observed phys- 
ical characteristics of the main phases we are tracing: the warm 
diffuse gas, which is widely distributed with low density, and 
the cold atomic phase which is composed of smaller and denser 
clouds. 



3.4. Data adjustment. 

The extinction data we used are not complete, in the sense that 
the more extinguished stars are not observed because of magni- 
tude selection. In the core of dense clouds, extinction reaches 10 
mag (Cambresy, 1997) with typical cloud sizes lower than a few 
pc. Thus, we underestimate the opacity in the darkest regions of 
the sky and the star sample used in this study allows us to de- 
tect only diffuse clouds. Moreover, the real density fluctuations 
occur at different scales and the smoothing length allows us to 
detect only typical dust cloud sizes larger than 5 pc which cor- 
responds to the smallest correlation length used in Eqn. |29] The 



spatial fluctuations smaller than this smoothing length are not 
detectable. 

Some extinction data are not compatible with our model (Fig|2 1, 
more particularly stars with large observed extinction. The co - 
ors of these stars are probably affected by structures smaller than 
the correlation length and cannot be taken into consideration in 
our model. 

For this reason, Figure|2]shows asymetric behaviour: some large 
observed extinctions are not correctly adjusted by our model and 
are underestimated. An other possibility is the presence of out- 
liers due to the fact that peculiar calibration stars could not be 
detected during selection step. 

4. Results 

The result of the inversion procedure is a 101^ data cube, cor- 
responding to a 500x500x500 pc^ cube in physical space, each 
individual 5x5x5 pc^ pixel being attributed a volume opacity. In 
this section we present the main characteristics of this opacity 
distribution, along with the results of a similar inversion of 1670 
interstellar neutral sodium (Nal) column densities measured 
towards stars within 300 parsecs. The sodium data are those 
presented in Welsh et al, (2009). At variance with the inversion 
shown in Welsh et al, performed for a unique kernel of 20 
pc correlation length, the present sodium inversion has been 
made based on the same data but for exactly the same inversion 
parameters as for the opacities, i.e. two kernels with correlation 
lengths of 5 and 15 pc respectively. One may expect some small 
differences from the Welsh et al maps, especially at the smallest 
scale, since the present inversion allows for smaller clouds. The 
choice of identical kernels is motivated by the need to eliminate 
one of the source of potential differences between sodium and 
dust opacity. 

This parallel analysis has two main reasons. First, both neutral 
sodium and opacity probe dense interstellar clouds and it 
is interesting to compare the distributions from the point of 
view of the ISM itself. The second reason, more important 
at this stage of our inversion studies, is technical. It is very 
informative to compare two inversions based on two completely 
independent data sets, because it provides some tests of the 
inversion, especially about the choice of the kernels in relation 
with the target distribution and its statistical significance. As 
we discussed earlier, the set of target stars for the opacity 
study is characterised by a non homogeneous distribution, with 
some regions particularly well covered, such as the Southern 
hemisphere region at 1< 20° and 1> 200°. Also, the density of 
stars decreases strongly beyond a distance varying between 100 
and 200 pc depending on the directions, due to selection effects 
linked to very strong extinction. The set of sodium target stars 
on the other hand is more regularly distributed with distance, 
but is a factor of four smaller. It is interesting to see the effect of 
such different grids on the large scale distribution of the inverted 
quantity. 



4.1. Integrated opacities 

4.1 .1 . Comparison with Sclilegel et al. full-sky dust maps 

An interesting comparison can be done using the reddening 
maps of Schlegel et al. (1998). These two maps (for the Northern 
and the Southern galactic hemispheres) give the total reddening 
E(B - V) deduced from dust emission measurements and a tem- 
perature correction. Using our extinction cloud model, we are 
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Fig. 6. Integrated opacity (resp. Nal column density) along 250 
parsec long lines-of-sight based on the inverted cubes, here rep- 
resented in Aitoff coordinates. 

able to estimate, at high galactic latitude, the total expected ex- 
tinction, integrating Eqn. |8] from to infinity (in practice, the 
integration is stopped at 250 pc). In order to compare directly 
the results of our model with those of Shlegel, we use the same 
grid representation and interpolate the integrated extinction on a 
similar sampling. 

Figures |3] and |4] show the comparisons separately for the north- 
ern and southern galactic hemisphere. This shows clearly that 
the extinction results are compatible: the extinction structures 
recovered from stars are underestimated according to the fact 
that cloud cores are not tracked by our sample. 

4.1.2. Comparison with integrated Nal columns 

A first global comparison between the two inversions can be 
done by integrating back the cube values, i.e. the volume opaci- 
ties and Nal densities, from the center to the surface of the cubes. 
Starting with the inverted cubes we have computed, for a regu- 
lar grid of directions and for all the sky, the total opacity and 
the total Nal column integrated from the Sun up to a distance of 
250 pc. The results are displayed in Fig|6]for a grid of directions 
evenly spaced by 2.5 deg. 

A first look at those figures reveals immediately the local galac- 
tic disk "warp", clearly visible here in both the dust and gas dis- 
tributions, with a general trend of concentration above the plane 
towards galactic longitudes 325 and 75 degrees, and below the 
plane elsewhere. This corresponds to the well-known Gould belt 
distribution of nearby O-B stars and gas. Such a clearly recon- 
structed "warp" in the two maps is a very encouraging result 
since the "a priori" solutions in both cases are plane-parallel and 
symmetric with respect to the plane. The similarities between 
the two integrated quantities are also obvious when consider- 
ing the large scales. This demonstrates that the stellar grids are 



Spatial distribution of interstellar dust in the Sun vicinity 

dense enough to reveal the main concentrations of dust and gas. 
A closer inspection reveals some differences at the smaller angu- 
lar scales, especially in the very dense regions. Such differences 
will be more clearly visible in the 2D cuts through the cubes 
presented below, and may correspond to (i) inaccuracies of the 
inversion due to a too coarse stellar grid, (ii) actual differences 
between gas and dust, (iii) inaccurate data. The detailed study of 
all these discrepancies is a very long task which is progressing 
along with the inversion of increasing amounts of targets. 
A clearly revealed feature is the marked difference between the 
integrated opacity at the North and South poles. This is to re- 
late to the North-South asymmetry already found by Schlegel et 
al (1998) in their reconstructed dust extinction maps based on 
COBE/DIRBE and IRAS data. That there is significant dust and 
a high dust/gas ratio in the North has also been inferred from ex- 
tinction and HI data by several works (Knude and H0g, 1999). 
Interestingly, the North-South asymmetry is not obviously 
present in the case of the sodium columns, and the comparison 
between the two maps suggests that the dust/gas ratio is differ- 
ent in the two polar areas within 250 parsecs. Such a difference, 
which should have some consequences in the context of dust 
emission models, certainly deserves further studies. 

4.2. Cuts in the 3D opacity and density cubes, identification 
of structures 

In this section we show the volume opacity and Nal density 
within several planes, i.e. different cuts in the 3D inverted cubes. 
For each map is shown the continuous value of the volume opac- 
ity (resp. density) with some iso- volume opacity (iso-volume 
density) contours drawn to help visualise the dense areas and the 
voids. Also shown are the stars which are close to the plane and 
have contributed to the adjustment. The criterion for keeping the 
stars is identical to the one presented in Sfeir et al (1999). The 
colour scale for the stars corresponds to the measured integrated 
opacity (resp. column density) between the Sun and the star. This 
allows to figure out how many stars have constrained the inver- 
sion and which stars have most influenced the local opacity (den- 
sity). 

Also shown in the maps are some markers of the reliability of 
the inversion process, which is directly linked to the number of 
targets and their spatial distributions: areas with averaging in- 
dex lower than 0.5 are marked by black dots. In almost all maps, 
the inversion is well constrained within 100-150 parsecs, while 
beyond this distances there are large differences according to 
the longitude or latitude. The most constrained plane for neu- 
tral sodium is the Gould belt plane, since there are numerous 
nearby early-type stars along the belt, which are ideal targets. 
High-latitude areas are not well sampled, at the exception of the 
poles themselves or well known "windows" of low columns. 

4.2.1. Vertical planes 

Figures|7]to[T0]display the inverted volume opacity and Na den- 
sity in four vertical planes containing the Sun at 45deg from 
each other, starting with the meridian plane, i.e. the vertical 
plane which contains the galactic centre and anti-centre direc- 
tions. Such vertical planes are the most appropriate cuts in the 
3D cubes for revealing the limits of the inversion process linked 
to the data sets, since the space density of the targets decreases 
strongly with the distance to the plane. Indeed, the influence of 
the "a priori distribution" can be easily discerned in those maps 
in all unconstrained areas, under the form of a smooth decrease 
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Fig. 7. Volume opacity (left) and Nal volume density (right) within the meridian plane. Added are several iso-contours, as well as 
the locations (triangles) of the target stars close to the plane (see text). The marker colour represents the integrated opacity (resp. 
column) towards the star. Small dots are plotted in areas where there are no constraints on the inversion, and the volume opacity 
(resp. density is equal to the a priori value. 
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Fig. 8. Same as figure [t] within a vertical plane at 45deg from the meridian. 



of the opacity (density) as a function of the distance to the plane 
and horizontal iso-densities. This is a visible criterion of lack of 
constraints. 

On the other hand, the juxtaposition of the same maps drawn 
from opacity and neutral sodium also clearly reveals the limita- 
tions of the 3D reconstruction. In some locations there are some 
differences between the centroids of the dense areas, or the ab- 
sence of a feature in one of the two maps. Those shortcomings 
are related to the target scarcity in some areas, or may also be 
related to parallax uncertainties. 

Despite those limitations, the Local Cavity is clearly defined in 
each of those vertical planes, and as noted in our previous works, 
is significantly inclined with respect to the plane. The global in- 
clination reaches maxima in the 0-180deg meridian plane and in 
the 315-135deg vertical plane, and the maps suggest a compli- 
cated structure with a number of low-density chimneys linking 
the cavity to the halo. 

Some nearby cavities are revealed by those vertical maps, almost 



in all cases in both extinction and gas maps when there are con- 
straints in those areas, but frequently the dust or gas distributions 
around the cavities differ slightly or more significantly, due to the 
target limitations. Work is in progress to identify these cavities 
in radio and X-ray maps and will be the subject of a foithcoming 
paper (Raimond et al, 2010). We have identified a few structures 
in the meridian plane in Fig|7] The distances to Taurus and pOph 
clouds will be discussed in the next section. The nearby Corona 
Australis concentration of clouds at -20deg is present but there 
are some substantial differences in the structure derived from ex- 
tinction and the sodium inversion. These differences may be due 
to the scarcity of targets stars in this area, and the small size of 
the CrA molecular region. The distance to the densest area in 
this direction is 135 pc from extinction, and about 165pc on the 
sodium map. The former distance corresponds to previous esti- 
mates, and the latter to the estimate by Knude and Hog (1998). 
More data are required in this area. 
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Fig. 9. Same as figurel?] within the rotation plane. 
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Fig. 10. Same as figurel?] within a vertical plane at 45deg from the rotation plane 



4.2.2. The Gould belt/Li ndblad ring plane 

In opposition to the vertical planes, the galactic plane and the 
Gould belt/Lindblad ring (GB) plane (see Poppel (1997) for a 
recent review on this structure) are among the most favourable 
regions for the inversion from the point of view of target 
distribution. For the extinction stellar dataset, the galactic plane 
is the most appropriate, but the Gould belt plane, being only 
slightly inclined, also contains a high density of targets. For the 
sodium absorption data set, it is the most appropriate, since it 
contains a large amount of bright, early-type stars that constitute 
ideal targets. On the other hand, the GB plane is highly struc- 
tured since it contains most of the dark clouds, star-forming 
regions and HIl regions. Although this complexity makes the 
inversion more difficult, it is interesting to attempt a mapping 
within this plane. We have thus chosen to display the inver- 
sion results within this particular plane, and the corresponding 



cuts in the two inverted cubes are displayed in Figures 1 1 and 12 



The Gould belt/Lindblad ring structure has been the subjects 
of extensive studies, and we have chosen here for the definition 
of the plane the coordinates calculated by Perrot and Grenier 



(2003), based on the distribution and kinematics of OB associ- 
ations and HI and H2 clouds. It is generally admitted that the 
belt, whose origin is still controversial, is about 600 pc wide 
and that the Sun lies within about 150 pc from its closest part, 
the Scorpius-Ophuchus complex of OB associations and dense 
clouds. The opacity and gas maps both reveal this complex 
area, and locates the East and West Ophiuchus clouds and 
the Lupus clouds from 100 to 165pc, in agreement with the 
recent work of Snow et al (2008) from absorption data (122 
±8 pc for the densest part of pOph and 133-150 pc for its 
western part), and Lombardi et al (2008) from reddening data 
(119 ±6 pc for the densest part of pOph and 155 ±8 pc for 
Lupus). The advantage of the present map compared to these 
two specific studies is that it provides a global view of all this 
region and the continuity between these regions. Similarly the 
map reveals the structure of the Coalsack (its southern part only 
is within this plane) and Chameleon dark clouds at distances 
similar to the previous results of Knude and Hog (1998). These 
authors found between 100 and 150 pc for South Coalsack, 
which corresponds well to the elongated dense area in our map, 
and two separate dense regions for Chameleon, the densest 
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Fig. 11. Volume opacity within the Gould Belt plane (orientation taken from Perrot and Grenier, 2003 



at about 150 pc and the second much closer, starting at about 
80 pc. This again corresponds well to the double structure 
we find in this direction. Opposite to this dense region, is the 
Taurus dark cloud complex, at about the same distance (150pc), 
in agreement with previous work on the young objects and 
molecular clouds in this area ( [Bertout et al.(1999 ) ). The star 
distribution allows to detect beyond Taurus the closest part of 
the giant Orion-Eridanus superbubble, along with some nearby 
condensations at its periphery. The Local cavity is opening to 
another giant cavity, the GS238+00+09 superbubble (Heiles, 
1998), adjacent to Orion-Eridanus. Interestingly, the Local 
Cavity has two extensions at l=85deg and 1=240-250 deg in 
this plane which are nearly perpendicular to the axis defined by 
the densest regions Taurus and Sco-Oph. The longest extension 
coincides with the well-know Canis Major extended void. 

Again, the comparison between the two distributions shows 
a rather good agreement about the global structure, which shows 
that the two datasets are good enough for such a large-scale 
determination. On the other hand, there are differences in the 
detailed structure. Work is in progress to identify the molecular 
complexes detected in radio data, and preliminary results 
confirm that, while some of the smallest isolated clouds are 
missed due to insufficient sky coverage, the main concentrations 
are retrieved. Further investigations are necessary in those cases, 
both the opacity and sodium star distributions are generally 



not sufficient to answer The most important diff'erence is at 
l=235deg, at a distance of about 200 pc. Here the dense con- 
centration detected in sodium does not seem to have a similarly 
strong counterpart in opacity. It is the contrary in the opposite 
direction, at about 120 pc. Is this due to a strongly variable dust 
to gas ratio, or to an artefact? 

From those maps and the distribution of the quality index, 
it is clear that there is a lack of targets beyond the densest con- 
centrations. E.g. it is not clear where the dust and gas densities 
drop beyond the Sco-Oph and Lupus clouds, or beyond Taurus. 
This is an important aspect, since the precise localization of the 
nearby cavities, potential X-ray emitters, would be a valuable 
advance. 

5. Conclusions 

We have built a carefully and conservatively selected list of 
opacity measurements based on Stromgren photometry, re- 
stricted to nearby stars possessing a Hipparcos parallax and lo- 
cated within 300 pc. We have inverted this opacity database 
(6400 targets) together with the interstellar sodium absorption 
database (1650 targets) presented in Welsh et al (2009) by means 
of an inversion tool specifically adjusted to those data. We have 
compared the resulting 3D dust and gas distributions in different 
ways, including maps of integrated opacity and density as well 
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Fig. 12. Same as figure[TT| for the neutral sodium volume density. 



as planar cuts within the 3D cubes. 
The main conclusions of those inversions are: 
-Integrals over 250 parsecs of the reconstructed opacities and 
sodium densities show very similar patterns over the sky, and re- 
veal in a similar fashion the local warp of the plane, a warp link 
to the Gould belt-Lindblad ring global structure. 
-The comparison with the integrated opacities from Schlegel et 
al. (1998) shows that the present results are compatible with 
those obtained from dust emission maps. 

-The extinction maps suggest a higher dust content at high galac- 
tic latitude in the Northern hemisphere compared to the South, 
while the gas maps do not. There is thus evidence for a higher 
dust to gas ratio in the North, in agreement with previous evi- 
dences (Knude and Hog, 1999) 

-Both tracers clearly reveal the Local Cavity and its contours are 
similar, despite the totally independent, and indeed very differ- 
ent data sets. 

-In general the large scale structures are in good agreement ev- 
erywhere there are constraints from the data. 
-At smaller scale there are some differences, that may be real but 
in most cases are likely due to the limited target coverage. This 
is clearly seen in the maps, especially far from the plane in the 
vertical cuts. Work is in progress on those discrepancies. 
The results of those comparisons are encouraging with respect to 
future work by showing that such inverse methods globally pro- 
vide similar structures despite the total independence of the two 



data sets and the strong differences in the target distributions. 
They demonstrate that the production of maps, with a precision 
of, say, 5 parsecs is not unrealistic. The extensive parallax and 
opacities databases from the forthcoming ESA mission GAIA 
and follow-up observations should allow to obtain such a pre- 
cision in the solar neighbourhood and to extend such maps to 
much larger distances. 
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